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Abstract. We characterize the force state of shear-loaded granular matter by 
relating the macroscopic stress to statistical properties of the force network. 
The purely repulsive nature of the interaction between grains naturally provides 
an upper bound for the sustainable shear stress, which we analyze using an 
optimization procedure inspired by the so-called force network ensemble. We 
establish a relation between the maximum possible shear resistance and the 
friction coefficient between individual grains, and find that anisotropics of the 
contact network (or the fabric tensor) only have a subdominant effect. These 
results can be considered the hyperstatic limit of the force network ensemble and 
we discuss possible implications for real systems. Finally, we argue how force 
anisotropies can be related quantitatively to experimental measurements of the 
effective elastic constants. 



PACS numbers: 45.70.Cc, 05.40.-a, 46.65.+g 

1. Introduction 

An assembly of cohesionless granular matter, in which there is no attraction between 
grains, can only exist when held together by an external pressure pp. The distribution 
of these confining forces throughout the material is a complex process that involves a 
highly inhomogeneous network of contact forces 130 EE]- Force networks as shown 
in figure^ are typical for a broad variety of amorphous systems like foams, colloids 
and emulsions, and play a crucial role for understanding the macroscopic mechanical 
properties [HlEIIHl. 

A robust feature of these "jammed" materials is that they can sustain a certain 
amount of shear stress before failure jHl EH HH There are many aspects that 

influence this shear resistance or internal friction of a granular material. A well known 
example is that wet sand can sustain much larger shear stresses than dry sand, due 
to the presence of attractive liquid bridges ^1 . The strength of the assembly is 
also enhanced by increasing intergrain friction or roughness of the grains. However, 
if one slowly increases the applied shear stress and follows the evolution of contact 
forces and grain locations, one encounters very complex collective phenomena. Before 
the system yields as a whole, there are non-adiabatic precursor events such as local 
rearrangements due to instability of subsets of grains |15[ I16L [T7| . This will induce 
changes in fabric and coordination number, and it has remained a great challenge to 
understand how this couples back to the stress state JJl^lElEIl- The conventional 
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Figure 1. (a) Force network obtained from a numerical simulation of a 
strongly hyperstatic frictionless granular material (coordination number z = 
5.75) subjected to a shear stress r = crxy/crxx, using the "force network 
ensemble" 5 21 . The thicknesses of lines represent the strength of the forces, 
(b) Corresponding average contact force as a function of contact angle {j> for this 
ensemble of force networks. Increasing the shear stress yields a modulation of 
f{4>) that is accurately described by the form 1 + 2Ts'm2if> + 62 cos 40. For large 
stress, f{<j>) approaches the limiting curve predicted by 1191 . 

tool to study this problem is direct numerical simulation of the particle dynamics. 
While this provides valuable information on the micromechanics of sheared granular 
materials, it remains very difficult to distinguish the relative importance of contact 
and force anisotropics: both are evolving simultaneously. 

Recently, a different strategy based on "ensembles of force networks" has been 
proposed to address this problem |2] . In this approach one investigates the statistics 
of all possible force network configurations that arc mechanically stable, for a single 
fixed packing geometry of grains 1,5, 22 , ,23, . This allows explicit separation of the 
effects of forces and texture, e.g. by studying the force ensembles for packs of different 
fabric and coordination numbers. Consistent with direct simulations '161, the ensemble 
showed that packings close to the minimum isostatic coordination number can hardly 
support any stress whereas strongly hyperstatic packings (much higher coordination 
numbers) can sustain much more stress |21j . At the same time, the theory also provides 
a very realistic description of the force anisotropy, force fluctuations and response 
function EH El 123 El 123 ■ 

In this paper we derive upper bounds for the shear load of cohesionless granular 
media, for varying intergrain friction and fabric anisotropy. The analysis is based on 
an intriguing observation made in the force network ensemble for strongly hyperstatic 
packs. The most elementary manifestation of force anisotropy is the modulation of the 
average force f (0), as a function of the contact orientation 0. From figure^, obtained 
from a two-dimensional frictionless system, it is clear that the force anisotropy is 
limited by the requirement that the normal component /„((/>) > for all (f>. This is due 
to the repulsive nature of the contact forces which require all /y > 0. Indeed, it was 
found that this simple criterion provides a very good approximation of the maximum 
shear stress achieved in the force network ensemble, in the limit of hyperstatic packs 
|21|. It also gives a good prediction for the limiting form of f{<j>) upon approaching 
the maximum load. 

The quantity f(0) thus provides crucial information on the force anisotropy and 
has recently been accessed experimentally for the first time [5j. To generalize the 
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Figure 2. (a) Definition of contact orientation rpij. (b,c) Illustration of the 
(average) bias of normal and tangential forces onto a particle due to the imposed 
shear stress. 

arguments above to frictional systems, we have to translate the physical constraints 
for all individual contact forces, 

• normal forces are purely repulsive, i.e. /„ > 0, 

• tangential forces ft obey Coulomb's law of friction, \ft\ < fJ- fn-, 

to constraints for f(0). Here, /x is the Coulomb friction coefficient of the individual 
contacts. We show how this yields nontrivial predictions for the maximum stress by 
finding the extreme forms of i{(j>)- While these maxima can probably not be reached 
in real systems, they describe the strongly hyperstatic limit of the force network 
ensemble, making it a well-defined analytical tool to investigate the influence of the 
micromechanical parameters on the effective macroscopic friction. 

The optimization approach followed in the present paper is in the same spirit 
as the force network ensemble, in the sense that it deals with the question if a force 
network with given parameters can exist or not, rather than with the actual evolution 
of a sheared granular system. On the other hand, while force balance on each grain 
is explicitly taken into account in the force network ensemble, it plays no role in our 
analysis of f (0). This is allowed in the strongly hyperstatic limit, but not for nearly 
isostatic packs. In the latter case the force networks are to a large extent determined 
by the conditions of force balance [Sj I^Hl ■ 

The relevant macroscopic quantity is the deviatoric stress, defined as r = 
((Ti — 0'2)l{o\ + 172), where the Ui denote the principal values of the stress tensor. 
We work in the coordinate frame where Oxx = Cyy, so we can express the deviatoric 
stress as T = dxyjoxx- For cohesionless systems one may naively expect that the 
ultimate shear stress is achieved when one principal direction becomes tensile, e.g. 
(72 = 0, which would lead to r^ax = 1- By invoking realistic structures of f((/)) for 
granular packs, however, we show that the physics is in fact much more subtle, and 
that the real maximum is typically much lower than unity. Our main findings are that 
the effect of friction between grains is only mild: a typical Coulomb friction coefficient 
of /i = 0.5 increases T„iax by only 16% as compared to the frictionless case. We also 
find that realistic anisotropics in the contact fabric hardly increase Tmax and thus seem 
to play a subdominant role in real systems. 

The paper is organized as follows. We first define the quantities that will be used 
to analyse the bounds, defined from the microscopic structure of grains and contacts, 
in relation to the macroscopic shear stress in sectional The third section addresses 
the simplest case of isotropic, frictionless packings, and shows how the granularity of 
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the material affects tlie maximally supportable shear stress. In sections 01 and [3 we 
explore the effects of intergrain friction and fabric anisotropies on the maximum shear 
stress. We conclude in sectional where we argue how our approach can be applied to 
problems of anisotropic elasticity. 

2. From contact forces to macroscopic stress 

We consider two-dimensional packings of discs, so that the orientations of the contacts 
between particles i and j can be characterized by the angle (figure 

= arccos (g^) • (1) 

Here r^ denotes the position vector of particle i and Xi its x-coordinate. The key 
quantity that we will use to determine the maximum shear load is the average force 
f (0) carried by all contacts of orientation (j). Since bond directions have no polarity, 
the angle only assumes values < (f) < tt, and the period of f((^) is tt. One can relate 
f (0) to the stress tensor cr, using |2H| 

o-a/3 = ^(f„)a(ry)^ = ^ f„r^ . (2) 

{ij} 

Here a, f3 label coordinate axes, Nc is the number of contacts in the (two-dimensional) 
volume V, Vij = Yj — r^, and is the force exerted on particle j by particle i. The 
bar indicates average over all forces in V . We decompose the force vector in a normal 
component, fn,ij, and tangential components, /t,y , as 

(fii)^ = fnAj COS - ft^ij sin (f)ij (3) 

(fij)a = fn,ij sin4>ij + ft,ij cos . (4) 

The sign conventions are such that repulsive forces have positive fn,ij, while the 
tangential component is positive when pointing "counterclockwise" with respect to 
particle i (see figure 0). 

For large enough packings we can express the stress tensor in a statistical form, 
evaluating the average faT/j from the probability to find a contact with force fij 
and centre-to-centre vector r^j . In terms of normal and tangential components, and 
with the observation that the forces are uncorrelated to the interparticle distance 
\rij I 1^1211, this involves the joint probability P{fn, ft, 4>)- We can exphcitly factorize 
the contact angle probability ^{(p) to write 

p(/„,/t,0) = $(0)P4/„,/t). 

The distribution P^{fmft) has been introduced recently |21| and is properly 
normalized to unity. Hence, the probabilistic form of the stress tensor reads 



fN f r°° 

= — ^ / d(/.<i>(0) / d/„ / d/t 
^ Jo Jo J-oo 



xP0(/„,/t)[/nAA„0 + /tr„/5] 

m r 
y Jo 



fN - - 

d^^c^) [Um^p + ftWT^p], (5) 
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where r denotes the average interparticle distance. The projection factors have been 
cohected in tensors Afap and 7^/3, written in matrix notation as 

/ cos^ (/) cos (/> sin (/) \ , , 

^ \^ COS (p sm (j) sm <P J 

J. _ f — COS (j) sin (p — sin^ 4> \ f'j\ 
"'^ ~ cos^ (j) COS (/) sin y ■ ' 

In the remainder, we wiU set the prefactor fNc/V = 1 in 

The above analysis allows computing the stress from f((/)). To relate the force 
anisotropics to the shear stress, a common trick is to expand f(0) in a Fourier 
series [13 IHl Ei 

N N 

fn{4') — ^ flfc sin 2k<j) + ^ bk cos 2k(j) (8) 
k=l k=a 
N N 

fi{(f>) = ^ Cfe sin 2k(f> + ^ dk cos 2k(t) . (9) 

k=l k=Q 

Because we are working in the frame where a^x = fyy, the principal axes of stress 
point in the directions (1, 1) and (1, —1). These directions must then be lines of mirror 
symmetry, as is illustrated in figure |21 which makes that ak = dk = for even k and 
bk = Ck = for odd k. 

For the moment, until the last section of the paper, we will consider the case 
where the fabric is isotropic so ^{(p) — I/tt. In that case, inserting equations (|8I9|) in 
((nj and integrating yields 

axx = bo/2 (10) 
ayy = bo/2 (11) 
axy = {ai+di)/i . (12) 

All higher order terms in the expansion yield zero upon integration. Our main interest 
is the deviatoric stress, so we are free to choose the pressure scale as a^x = <^yy = V^i 
so that / = 6o = 1 and 

r^^. (13) 

This relation reveals how an applied shear stress can be sustained through anisotropics 
in both the normal and frictional forces, via ai and di respectively [22] • The strategy 
will be to explore the physical limitations of Oi and di, which will provide a bound on 

T. 

Note that due to the linearity of ((2Jl, the stress only couples to the first moment 
of the force distributions. This means that the stress does not depend on details of 
the probability -P(|f|) jH ^ HI [Z| |H| . Also, the stress tensor contains no information 
on "force chains" , i.e. the tendency of large forces the align in a correlated way. (|2Jl 
does not invoke products of different so that spatial force-force correlations do not 
come into play. 



3. Frictionless packings 



We start out with packings of frictionless particles, for which the problem of the 
maximum possible deviatoric stress is relatively straightforward. In this case a bound 
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on r emerges from the purely repulsive nature of the forces: all contacts have fn^ij > 0, 
so certainly the averages should obey 

W) > , (14) 

for all values of (p. This condition obviously forms a serious restriction on the amplitude 
of the force anisotropy, and we show how this yields an upper limit on r. We show that 
this maximum value Tmax depends on the number N of terms included in the Fourier 
series of (jSJ: even though the higher order terms do not contribute to the stress, they 
enable f^^ij to reach more extreme shapes. We first explore the case N — 2, which is 
the relevant case for real systems (see also |Appendix A| |. Then we discuss the problem 
for arbitraty N, which provides some additional insights. 



3.1. Realistic /n(0) 

For the frictionless case we have from 

ai = 2t, di = 0, (15) 
so that (|SJ) becomes 

/„(0) l + 2Tsin20 + 62Cos4(/)H (16) 

Let us first consider the simplest case, in which we truncate after the lowest anistropic 
term, i.e. fn{(f>) = 1 + 2Tsin20, so that r is the only free parameter. For positive r 
this function has a minimum in = 37r/4, which touches /n(</') = for 2r = 1. Hence, 
Tmax = 1/2. In the same way we could derive Tmin = —1/2. From now on we will only 
consider positive r, without loss of generality. 

From the numerical result of figure^ it is clear that the modulation does not 
stay symmetric around /n(0) = 1 for large values of r. This indicates a significant 
contribution of the type cos This is the highest order term that we can expect to 
be relevant in granular matter, because steric exclusion between the grains sets a lower 
limit to the width of peaks in /n(0)- This is explained in more detail in [Appendix A| 

We thus focus on the case = 2 in the expansion ^ . The optimization problem 
involves two free parameters, r and 62- We are free to vary &2 in such a way as to 
facilitate a maximum r, under the constraint that /n(<^) > for all (j). We furthermore 
demand that /n('/') evolves monotonically between the principal directions at (/> = 7r/4 
and (j) = 37r/4. This implies that the minimum of fniff) should stay at 4> = 37r/4, as it 
is the case in the numerics of figure Physically, this means that the first contacts 
to break are those oriented in the direction in which the material is stretched. While 
the expansion ensures that there is an extremum in </> = 37r/4, this extremum only 
remains a minimum if 

= 8r + I662 > . 
The value of /n(</') in this minimum must satisfy 




/n(37r/4) = l-2T-62>0 
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Figure 3. The optimized /n(0) for = 2 (solid), = 10 (dashed). The figure 
in the corner illustrates the corresponding average force exerted on a particle for 
various contact orientations (A'^ = 2). 



This defines a linear program with parameters r and 62, two inequahties, and the 
objective to maximize r. The solution is easily found to be | 
2 

Tmax = g (17) 
h = - ^ . (18) 

figure^ illustrates the relevance of this bound for strongly hyperstatic packings: the 
numerical fn{4') (taken from [21]) indeed approaches the limiting form (dashed curve) 

/n(0) = l + ^sin20-icos4<?!). (19) 

So indeed, the system is able to organize the forces in such a manner as to optimize 
the sustained shear stress. 



3.2. The limit N —^ 00 and tensile stresses 

Although Fourier terms with fc > 3 do not play a role in granular systems, it is 
insightful to study the general case where we truncate the series at arbitrary order N. 
The expansion now contains N free coefhcients that we need to fix, in order to optimize 
r. In the case = 2, we invoked the condition d^fn/dcj)^ > a.t cj) — 3tt/A to ensure 
that the minor principal axis remains a minimum of /n(0), and taking d'^fn/d(fP = 
turned out to maximize r. A nonzero fc = 3 term allows to also fix d^fn/dcj)'^ in such 
a way that r becomes even larger. The upshot of adding more terms is that we can 
make /n(0) as flat (and as close to zero) as possible around (p = 37r/4, where the 
contribution to the overlap with sin 20 is negative. Every time we have to verify that 

f Actually, because for these values of the parameters the second derivative vanishes in t/i = 37r/4, 
we have to check that the fourth derivative is positive to be sure that the extremum is a minimum, 
which is indeed the case. 
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the first nonzero Taylor coefRcient when expanding around (j) — Stt/A is positive, so 
that we are indeed dealing with a minimum. 

The general scheme is thus that adding the term of order k generates an additional 
condition that = 0. For general iV, this yields the following set of 

linear equations: 



& 



ai 



4,21 



for all I = 0,2,---,iV - 1 . (20) 



In [Appendix B| we show that this linear problem for the Fourier coefBcients can be 
inverted analytically, yielding a remarkably simple result for the maximum r, namely 

N 

w(iV) = . (21) 

Interestingly, this reveals an ultimate (mathematical) maximum shear stress 
Tmax = 1- We plotted the optimized /n(0) for various values of N (figure 0)): clearly, 
/n(0) evolves towards the extreme case of a Dirac (5-peak in the limit N ^ oo. As 
already mentioned in the introduction, this condition of t = 1 precisely corresponds 
to the point where the minor principal axis becomes tensile. We thus conclude that, 
due to the finite width of /n(0), the maximum stress for a frictionless packing lies 
well below the point where the global stress unavoidably develops a tensile direction. 
We have argued in [Appendix A| that this finite width is in fact due to steric exclusion 
between neighbouring grains. This illustrates how the discrete nature of the assembly 
has an important effect on global properties. 



4. Prictional packings: T„^ax{f^) 

The presence of frictional forces provides an additional degree of freedom to develop 
anisotropic stresses: Equation shows that the total deviatoric is the sum of 
the (lowest order) anistropies of normal and tangential forces. It is clear that this 
will enhance the ability to sustain a large external load. However, there is again a 
bound on the force anisotropics, now due to Coulomb's law for individual contacts, 
i-S- Ift.ijl < fJ'fn.ij, where /i is the microscopic Coulomb friction coefficient. Since this 
condition should hold for any pair of grains, it certainly holds for the averages: 

l/t(0)l < ■ (22) 

This condition is illustrated in figure 0] In this section we derive how Tmax depends 
on the value of /i, again using the Fourier expansions of /n(0) and ft{4') up to A'^ = 2. 

4.I. The optimization problem 
From H13|l we can express 

ai = 2t - rfi , (23) 

so that 

/„ (,/,) = 1 + (2t - di ) sin 20 + &2 cos 40 , (24) 
(0) = di cos 20 + C2 sin 40 . (25) 

Let us now take a value of r slightly above the frictionless limit 2/3. The prefactor 
in front of the sin 20 term can now be lowered due to di, i.e. due to the presence of 



Bounds on the shear load of cohesionless granular matter 



9 




Figure 4. The frictional optimization problem for /i = 0.28. The solid lines 
represent it /n (</>), the dash-dotted line /t((/))//i, which has to satisfy I22i . The 
figure in the corner illustrates the average force exerted on a frictional particle 
for various contact orientations. Note that these forces now have tangential 
components. 



friction. If we put r far above 2/3, however, one requires a relatively large di. The 
value of di is bounded by the condition of (|22|l . so that not all r can be reached. 

To determine the maximum value of t as a function of fi, we have to specify 
acceptable values of the higher order coefficients 62 and C2 , which do not contribute to 
the stress tensor. As we did in the frictionless case, we demand that the function /n(0) 
evolves monotonically between major and minor directions, so that it has only one 
maximum (in the major direction) and only one minimum (in the minor direction). A 
similar requirement is imposed on /t(0): the average tangential force only changes sign 
along the principal directions, see e.g. figure01 If this were not the case, the tangential 
forces would swap from clockwise to counterclockwise and back in between the major 
and minor directions. Such a spontaneous symmetry breaking would introduce a very 
nongeneric organization of forces within the packing. These conditions put bounds on 
the second order coefficients. 



2t - di 



<bo< 



2t- di 



(26) 
(27) 



4 - - - 4 
di ^ ^ di 

< C2 < . 

2 - 2 
numerically solved this 
of the parameters r, di 
equations I|22I26I27|I . The results are shown in figure Surprisingly, the dependence 
on /i is relatively weak. In the paragraph below we obtain the analytical result. 



We have 
possible values 



optimization problem by 
,62,22, within the ranges 



varying 
imposed 



all 

by 



1 + VT+V 



which is derived for 1 — < fu, < 1. However, figure El shows that this is a very 
good approximation outside this range as well. 
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Figure 5. The optimized t(^) for frictional packings with isotropic fabric. The 
circles are numerical data. The solid line is equation 1281 . The dashed line is 
the analytic result for < 1 — -^V^ ~ 0.29. The dotted line is the asymptote 
T (X 0.593fi as /I — + oo. The inset shows an enlargement of the small /i region. 

4-2. Analytical solution ofr^fi) 

The set of parameters t, di , 62 , C2 that corresponds to the maximum value of t 
obviously corresponds to functions ft{(f>) and ^/n(0) that are tangent in at least one 
point. We will now derive the optimal set of parameters for the case there is a tangent 
point in the interval < < tt, and then determine for what range of /i this is the 
case. 

In the interval tt], the cos 20 in /t(0) is positive and the sin 40 is negative. We 
want to have a c?i which is large (to facilitate larger r), and at the same time a ft{4>) 
which is as close to zero as possible (to stay away from violating (|22|l l. Therefore, the 
parameter C2 should have its maximum value C2 — rfi/2, as dictated by H27|) . Similarly, 
we want /n(0) to be as large as possible, and hence because cos 40 is positive in the 
considered interval, 62 should also be maximal. There are two upper bounds on 62, 
given by l|26|l and /n(0) > 0, the latter of which turns out to be the most restrictive 
one. This gives 

62 = 1 + di - 2r . 

Having eliminated 62 and C2, we have thus reduced the optimization problem to 2 
parameters, namely r and di. To find the parameters that maximize r, we take the 
equality sign in H22(l and demand that solutions are tangent points. For arbitrary 
values of the parameters, there is a tangent point in 0i = 37r/4 and two normal 
intersection points in 02 and 03. The intersection points should coincide to turn into 
a tangent point. Equating 02 = 03 gives a relation between the parameters which, 
after some lengthy but elementary algebra, can be written as 

Maximizing this expression for r(/Lt,di) with respect to di yields that the optimum is 




Bounds on the shear load of cohesionless granular matter 



11 



obtained for 



so that 



2^l 

di = 



1 + y/l + 

rmax(Ai) = ■ (28) 

Inserting these into the original equations gives the value of (j) where the curves touch, 

tan20=^^i±M^, (29) 
3/i 

which was demanded to be in the interval tt]. This requirement is met for values 
of n satisfying 

(0.29 w) 1- ^^2 < n<l . (30) 

This validity is indeed found numerically in figure[3 Note that this range of Coulomb 
coefficients is the most relevant for real granular materials |3()| . 

For smaller ^ the tangent point is below (p = where the cos 40 is negative, 
so 62 should now be as small as possible: 62 — — (2r — d)/A. The resulting quartic 
equations can be solved using computer algebra, yielding a lenghty expression (not 
shown), which is plotted as the dashed part of the curve in figureEl and which coincides 
with the numerical data. 

For > 1 the above analysis gives a tangent point in the interval [0,7r/12], so 
that the considerations that allowed us to fix 62 and C2 are no longer valid and we 
only have the numerical result of figure[5l When the Coulomb coefficient becomes very 
large, the sustainable shear stress is mostly due to the tangential forces. Because the 
right hand side of H22|l grows linearly with /i, the values of di and C2 can also scale as 
fj, when /i ^ 1. This leads to an to an asymptotic behaviour of Tniax(M) which is linear 
in /i. We numerically determined the parameters ai, 62, C2//i, di/ fi for the asymptotic 
/„((/)) and /t(0)/M by optimizing for di/fi instead of t. The resulting asymptote is 
r = 0.376 + 0.593^, which plotted as the dotted hue in figure 



5. Fabric anisotropy 

We now analyze the effect of a fabric anisotropy on the maximum shear stress. 
The structure of the contact network can be characterized using the fabric tensor 
■Pa/3 — TaT^/r^, and a fabric anisotropy shows up as a difference between the principal 
values of this tensor, Fi — F2. Analogous to the stress tensor, this difference is solely 
determined by the lowest order coefficient of the Fourier expansion of <&(0), the contact 
angle distribution, which we therefore approximate as § 

$(0) = -(1+Pisin2(/)) . (31) 

TT 

One easily shows that pi — 2{Fi — F2). The strategy is to again find the maximum 
possible value for r, but now for pi ^ 0. 

Let us note that this truncated form for $(0) is a good approximation for systems 
with a simple shear deformation history [291 131| , although more complex forms are 

§ The maximum shear stress is achieved when contact and force anisotropy are "in phase", so that 
we do not introduce a phase shift in the sin 2if) term of 13 li . 
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Figure 6. Maximum shear stress Tmax as a function of friction coefficient fi for 
pi = 0,0.3,0.5. 



encountered e.g. for packings formed under gravity As was demonstrated by 

Troadec et al. [23, the values for pi are bounded by the eflFect of steric exclusion 
between neighbouring contacts (in the same way as steric exclusion bounds the force 
anisotropy, see also [Appendix A\ . Indeed, numerical simulations show that pi < 0.3, 
while typically pi is of the order 0.1 j29| 1311. 

We thus insert the form l|31|l in the expression for the stress tensor yielding 

a^x = bo/2+piai/4 (32) 
ayy ^ ba/2 + piai/A (33) 

<Jxv^]{ai+di) + ^{2bo + c2-b2) . (34) 
4 o 

This time the second order coefficients do affect the value of the deviatoric stress 
directly. The idea of the optimization procedure, however, remains the same as 
in the isotropic case. In this section we present the results of the optimization. 
The frictionless case, which has an analytic solution, is discussed in more detail in 
[Appendix C| 

We again eliminate bo by setting axx = 1/2 in H32() . From (|34() we then have 

^ _ ai(2 - pI) + 2di +pi(2 - + C2) ^^^^ 

Working out the optimization for the frictionless case (c2 — di = 0), we get the 
analytic result 

_ 7pi+8 

^'""'^ ~ 8pi + 12 ' ^^^^ 
as is shown in [Appendix C[ The inclusion of fabric anisotropy hence leads to a small 
increase of the maximum deviatoric stress; for e.g. pi = 0.3 we get a 5% increase. See 
also the numerical data points on the vertical axis of figure jS] 

For the complete case, with anisotropic fabric and finite friction coefficient, we 
have used the same numerical optimization procedure as in section [4. II The result is 



Bounds on the shear load of cohesionless granular matter 



13 



shown for pi = 0.3,0.5 in figure El From this figure one can see that the combined 
effect of friction and anisotropy can roughly be seen as an addition of their individual 
effects. 

6. Discussion 

The repulsive-only interaction between dry grains causes the shear stress sustained 
by granular materials to be bounded. We have derived these upper bounds, Tmax, by 
finding the extreme shapes of the angle- resolved average force, f (0), that are consistent 
with Coulomb's friction law for individual grains. In the spirit of the Force Network 
Ensemble, the key question addressed in the present paper is "can such a force network 
exist?" , and we have argued that our analysis represents the strongly hyperstatic limit 
of the ensemble. 

The yielding process in real systems is much more involved and cannot be fully 
described by the existence criteria used in the paper. However, the analysis does allow 
to investigate the influence of the micromechanical parameters on the maximum shear 
stress, in such a way that different aspects can be disentangled. This provides useful 
complementary information to simulation methods, for which the packing texture 
cannot be controlled during loading. For example, we have found that an anisotropic 
fabric is not needed to sustain a large shear stress, and in fact adding fabric anisotropy 
hardly increases the maximum stress. This suggests that shear-induced textures 
observed in numerical simulations |29[ 13 Ij play a relatively passive role in the stress 
balance. 

Secondly, we analyzed the effect of the microscopic friction coefficient fi, and 
found that T,„ax does not increase rapidly with fj, in the physically relevant regime 
(figure|3J). One has to be very careful, however, to interpret this result. The presence 
of friction does have an important influence on the coordination number, which in 
turn affects the shear resistance: our analysis breaks down close to the isostatic limit, 
at which the force balance equations become dominant and Tmax drops to zero j21j . 
Consistent with our analysis, the internal friction measured in simulations of quasi- 
static shear flows displays only a mild increase is observed when fi > 0.4. On 
the other hand, t was found to decrease rapidly as 0, which we believe is because 
the system approaches the isostatic limit ||. 

Our main conclusion is thus that strongly hyperstatic packings of grains can 
support an amount of shear stress that is not very sensitive to friction or fabric 
anisotropy. For systems close to the isostatic limit, however, the dependence of the 
coordination number on any of these parameters has a more dramatic effect. 

Another interesting perspective is that one can use an expansion of f (</<), as 
presented in this paper, to estimate the effective elastic moduli of the system, denoted 
by Ca/s-yS- It has been shown experimentally that these become highly anisotropic 
when a system is pre-sheared |S1 EHl ■ The response to a localized vertical force on the 
free surface of a granular bed was found not to propagate along the vertical, but along 

II It should be noted that it is not obvious to determine how far from isostaticity (marginal stability) 
a system is. Two-dimensional systems of frictionless disks in principle become isostatic at ^ = 3, 
but when the microscopic friction coefficient ^ is small, many contacts are fully mobilized (i.e. the 
frictional force has its maximum value allowed by the Coulomb criterion) )35| . These contacts 
contribute less to the number of degrees of freedom in the force network, since the tangential 
force component is fixed by the normal force component [Tl 1351 . Therefore, these systems are less 
hyperstatic than one would think on the basis of the coordination number. 
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a direction that is tilted towards the major stress axis. This indicates a stiffening of 
contacts along the major axis that are responsible for the anisotropic elasticity. By 
fitting the experimentally measured stress profiles Atman et al. |36j obtained the ratio 
C1111/C2222 ~ 0.67 of the Young's moduli in the minor (1) and major (2) directions. 
This finding was remarkably insensitive to frictional properties and roughness of the 
grains. 

The stiffening of contacts can be explained through the nonlinear interaction 

between particles, which for frictionless Hertzian contacts |37) in three dimensions 
1/3 

is ^ij ^ 

f^J , where kij is the effective spring-constant of the contact. Assuming 
that the displacements of the particles are afhne — this is reasonable for the strongly 
hyperstatic packings considered in this paper — one can estimate the elastic 
moduh as PnUm 1 



This allows to systematically explore the effect of stress induced anisotropy in forces or 
fabric, by again using the Fourier expansions of /(</)) and $(0). For the simplest case of 
isotropic frictionless contacts, one already finds that C1111/C2222 = 1 — 8t/9 + O(t^), 
so that the experimentally observed ratio of about 0.67 is easily achieved for realistic 
values of r. 
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Appendix A. Steric exclusion and the coordination number 

We show in section |21 that optimizing r^jax is equivalent to making the maximum 
of f (0) around the major principal axis as sharp as possible. It has been reported 
from numerical simulations of realistic systems that in practice this maximum is not 
very sharp and that typically the highest Fourier contribution to f((/)) is of order 
cos4(/) |21[ I29j (see also figure H^)- This can be related to the coordination number 
as follows. Suppose we wish to maximize the forces around the major principal axis, 
(j) = 7r/4. On a typical grain this can be achieved by putting all stress on the contact 
closest to this angle. For a randomly oriented grain with 2 = 6, the closest contact 
lies roughly in the range of ±7r/6 around the principal direction. Let us imagine that 
it were possible for a grain to have z ^ 6. Now, the contacts closest to the major 
direction are spread within a much smaller range ±7r/z, so that much narrower widths 
of f(0) can be achieved. Of course, the situations described above are very extreme 
because in reality all contacts will be involved in the force equilibrium. In any case, 

^ The use of (f)^^^ instead of f^^'^ introduces only a small error since these functions differ mainly 
around the minor principal axis where the integrand is small anyway. 




(37) 
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these simple examples illustrate that the coordination number introduces a natural 
scale for the width due to steric exclusion effects. 

If the polydispersity is small, only very few grains will have more than six 
neighbours, and the peaks and valleys in {{4>) will therefore have a minimum width 
of approximately 7r/3. This is why a description using only a few terms in a Fourier 
expansion works, as was already clear from existing numerical results. 



Appendix B. Prictionless case at arbitrary order 

In this appendix we analytically solve the linear problem of (|20() . Due to the symmetry 
of /n(0) we can write 

/„(</>) = 1 + E* sin (^2fc0+^^il^^ ■ 



Comparing to lfTB|) . we find qk = (-1)'*^ ^^^^flfc for odd fc, while qk — (— l)'^^^6fc for 
even k. In particular, qi = 2r. The even derivatives can be written as 

^ ^ (-1)' J2 '?fc(2^)'' (^fc'^ + , (B.2) 

and since all sine terms evaluated at (f) — 3tt/A yield —1, we find for / 7^ 



^21 N 



4,21 



= (-l)'+i^g,(2fcf . (B.3) 

i-K 

— k=l 

The linear problem of H20|l can now be written in the form of a Vandermonde matrix. 



/ 1 



X2 



1 

2 

N 



X 



92 
<73 



V 



N-1 



( 1 \ 






(B.4) 



j j V / 



with Xk = 4fc2. The inverse Aj^ of this matrix can be expressed explicitly in terms of 
Lagrange polynomials as 01] 



N 

n 



Xj X}^ 



N 



(B.5) 



k=l,k^j ■' n=l 

We are interested in the solution for qi ~ 2t, which for l)B.4(l simply reads 



N 



Xk - Xi 



fc=2 



N 

n 

k=2 



2N 



(fc + l)(fc-l) 7V + 1 ' 



(B.6) 
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so that r = N/{N + 1), (PT|l . Similarly, the other follow from 
TT 



Appendix C. Details for the anisotropic frictionless optimization 

We start from the expression for the deviatoric stress in terms of all parameters and 
Fourier coefficients, equation (|35|) . where we put C2 — di — because the tangential 
force components are zero. This reads 

_^ ai(2~pj)+Pi(2-62) ^^^^ 

The physical constraints on ai and 62 are the same as for the isotropic case discussed in 
section IXn but the optimization target is now given by HC.1|) . instead of the r = ai/2 
we had before. The constraints follow from /n(0) and its second derivative being 
nonnegative at = 3tt/A: 

62 < l-ai(l+Pi/2) (C.2) 

62 > -ai/4 . (C.3) 

The solution to this linear program is then found to be 

ai - 4/(2pi + 3) (C.4) 

6o-3/(2pi + 3) (C.5) 

b2= - l/(2pi + 3) , (C.6) 



which corresponds to 

7pi 



8pi + 12 ' 

which is the result stated in equation (|36|l . 



(C.7) 
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